test_stockassessment

Setup

library(tidyverse)
#library(fishmethods)

Fmsy Equation

Test the Fmsy equation: Fmsy = r * (1 / (1+p)) - from the catchmsy function in R

fmsy = 0.2 * (1/(1+0.2))

Goals and Model

=======

#{.tabset}

##Setup

library(tidyverse)

##Fmsy Equation

Test the Fmsy equation: Fmsy = r * (1 / (1+p)) - from the catchmsy function in R

fmsy = 0.2 * (1/(1+0.2))

##Goals and Model

>>>>>>> e39bf25c59cba56a453798669eb224ec741eec47

Initial Stock Assessment Test Goal: Simulate error around estimates of intial biomass

Write the function to include error in the biomass estimates. In this experiment the b_err term represents the estimated biomass with error. The biomass is randomly selected from a distribution where the mean is the actual biomass from the model and the standard deviation is 5% of the actual biomass from the model.

The basic function that we’ll use in the rest of these tests:

#Write the function
sim_stockerror <- function(b, r, r_s, error, p, k, years){
  
  results <- data.frame(
    b = rep(NA, years), c = rep(NA, years), 
    year = 1:years, r = rep(NA, years), b_err = rep(NA, years)) #Setup the results dataframe 
  
  #Set the initial result for the outputs in year 1
   
  results$b[1] = b
  results$b_err[1] = rnorm(1, mean = results$b[1], sd = (error*results$b[1])) #simulate error in assessment
  results$c[1] = (r*(1/(1+p))) * results$b_err[1] #using f=Fmsy and b=estimate from initial stock assessment
  results$r[1] = r
  
  #Loop the model over the specified number of years
  for (t in 2:years) {
    
    results$b_err[t] = rnorm(1, mean = results$b[t-1], sd = (error*results$b[t-1]))
    results$c[t] = (r*(1/(1+p))) * results$b_err[t]
    results$b[t] = results$b[t-1] + (results$r[t-1] / p)*results$b[t-1]*(1 - ((results$b[t-1]/k) ^ p))-results$c[t]
    results$r[t] = results$r[1] * (1 + (r_s*(t-1))) 
    
   } 
  
  return(results)
}

#Test the results
results = sim_stockerror(b = 1000, r = 0.2, error = 0.05, p = 0.2,
                      years = 100, r_s = 0, k = 10000) 


plot(results$year, results$b)

<<<<<<< HEAD

plot(results$year, results$b_err)

Goals Tests 1-4

Apply the new function to a set of experiments that vary actual initial biomass. In these experiments catch is responding to underlying changes in actual biomass but with some error. - Each test will change only the standard deviation associated with the error term

Test 1

Standard deviation = 0.05

=======

plot(results$year, results$b_err)

##Goals Tests 1-4 Apply the new function to a set of experiments that vary actual initial biomass. In these experiments catch is responding to underlying changes in actual biomass but with some error. - Each test will change only the standard deviation associated with the error term

##Test 1 Standard deviation = 0.05

>>>>>>> e39bf25c59cba56a453798669eb224ec741eec47

n_experiment = 10 

#Input variables

experiment <- list(
  b = seq(1000, 10000, 1000),
  r = rep(0.2, n_experiment),
  r_s = rep(0, n_experiment),
  error = rep(0.05, n_experiment),
  p = rep(0.2, n_experiment),
  k = rep(10000, n_experiment),
  years = rep(100, n_experiment))

#Now apply to the function
results_test1 <- pmap(experiment, sim_stockerror)

#Plot results to see variation:
for(i in 1:n_experiment){
  plot(results_test1[[i]]$b)
  plot(results_test1[[i]]$b_err)
}

<<<<<<< HEAD

Results: Looks like regardless of initial biomass, with a error of only 5% biomass stays around Bmsy for all 100 years

Test 2

Standard deviation to 0.10

=======

Results: Looks like regardless of initial biomass, with a error of only 5% biomass stays around Bmsy for all 100 years

##Test 2 Standard deviation to 0.10

>>>>>>> e39bf25c59cba56a453798669eb224ec741eec47

n_experiment = 10 

#Input variables

experiment <- list(
  b = seq(1000, 10000, 1000),
  r = rep(0.2, n_experiment),
  r_s = rep(0, n_experiment),
  error = rep(0.1, n_experiment),
  p = rep(0.2, n_experiment),
  k = rep(10000, n_experiment),
  years = rep(100, n_experiment))

#Now apply to the function
results_test2 <- pmap(experiment, sim_stockerror)

#Plot results to see variation:
for(i in 1:n_experiment){
  plot(results_test2[[i]]$b)
  plot(results_test2[[i]]$b_err)
}

<<<<<<< HEAD

Results: About the same, regardless of initial biomass, with a 10% error biomass stays around Bmsy for all 100 years

Test 3

Standard deviation to 0.15

=======

Results: About the same, regardless of initial biomass, with a 10% error biomass stays around Bmsy for all 100 years

##Test 3 Standard deviation to 0.15

>>>>>>> e39bf25c59cba56a453798669eb224ec741eec47

n_experiment = 10 

#Input variables

experiment <- list(
  b = seq(1000, 10000, 1000),
  r = rep(0.2, n_experiment),
  r_s = rep(0, n_experiment),
  error = rep(0.15, n_experiment),
  p = rep(0.2, n_experiment),
  k = rep(10000, n_experiment),
  years = rep(100, n_experiment))

#Now apply to the function
results_test3 <- pmap(experiment, sim_stockerror)

#Plot results to see variation:
for(i in 1:n_experiment){
  plot(results_test3[[i]]$b)
  plot(results_test3[[i]]$b_err)
}

<<<<<<< HEAD

Results: About the same, regardless of initial biomass, with a 15% error biomass stays around Bmsy for all 100 years

Test 4

Standard deviation to 0.2

=======

Results: About the same, regardless of initial biomass, with a 15% error biomass stays around Bmsy for all 100 years

##Test 4 Standard deviation to 0.2

>>>>>>> e39bf25c59cba56a453798669eb224ec741eec47

n_experiment = 10 

#Input variables

experiment <- list(
  b = seq(1000, 10000, 1000),
  r = rep(0.2, n_experiment),
  r_s = rep(0, n_experiment),
  error = rep(0.2, n_experiment),
  p = rep(0.2, n_experiment),
  k = rep(10000, n_experiment),
  years = rep(100, n_experiment))

#Now apply to the function
results_test4 <- pmap(experiment, sim_stockerror)

#Plot results to see variation:
for(i in 1:n_experiment){
  plot(results_test4[[i]]$b)
  plot(results_test4[[i]]$b_err)
}

<<<<<<< HEAD

Results: About the same, regardless of initial biomass, with a 20% error biomass stays around Bmsy for all 100 years

Goals Tests 5-7

Use a wider range of error terms and systematically add in other variations: biomass, r, and r_s

Test 5

Varying biomass with a wider range of error terms

=======

Results: About the same, regardless of initial biomass, with a 20% error biomass stays around Bmsy for all 100 years

##Goals Tests 5-7 Use a wider range of error terms and systematically add in other variations: biomass, r, and r_s

##Test 5 Varying biomass with a wider range of error terms

>>>>>>> e39bf25c59cba56a453798669eb224ec741eec47

First create the correct list of error terms to input in the model:

#First try creating a dataframe with the sequences we want to repeat:
input <- seq(0.1, 0.55, 0.05)

#Create a list for the output
error_input = data.frame(error = rep(NA, 100))

#Loop the repetition of each input number
for(i in 1:100){
  if(i <= 10){
   error_input$error[i] = input[1] 
  } 
  if(i > 10){
    error_input$error[i] = input[2]
  }
  if(i > 20){
    error_input$error[i] = input[3]
  }
  if(i>30){
    error_input$error[i] = input[4]
  }
  if(i>40){
    error_input$error[i] = input[5]
  }
  if(i>50){
    error_input$error[i] = input[6]
  }
  if(i>60){
    error_input$error[i] = input[7]
  }
  if(i>70){
    error_input$error[i] = input[8]
  }
  if(i>80){
    error_input$error[i] = input[9]
  }
  if(i>90){
    error_input$error[i] = input[10]
  }
}

#This works but it seems very inefficient - is there a better/faster way to get this result?

Run the experiment:

n_experiment = 100 

#Input variables
experiment <- list(
  b = rep(seq(1000, 10000, 1000), 10),
  r = rep(0.2, n_experiment),
  r_s = rep(0, n_experiment),
  error = error_input$error,
  p = rep(0.2, n_experiment),
  k = rep(10000, n_experiment),
  years = rep(100, n_experiment))

#Now apply to the function
results_test5 <- pmap(experiment, sim_stockerror)

<<<<<<< HEAD